Estadística de Áreas y Autocorrelación Espacial

Geoestadística: Puntos en el espacio y variogramas para construir modelos :: Estadística de Áreas: regiones y vecinos para contrastar hipótesis de dependencia y para construir modelos.

Contenido:

  1. Vecinos, vecindades y pesos. Hoy
  2. Medidas de autocorrelación espacial. Hoy
  3. Modelos para medidas dependientes en áreas. La otra semana

Datos:

281 segmentos censales en EEUU, para 8 condados del estado de NY.

Vecinos en el espacio

¿Cómo se representa en forma matricial quién es vecino de quién? Paquetes disponibles: spdep, ade4 (nb2beig, neig2nb),

## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.1.1, released 2020/06/22
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/marce/Dropbox/UCR_2020/SP1649-Espacial/0.clases/data/NY8_utm18.shp", layer: "NY8_utm18"
## with 281 features
## It has 17 fields
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/marce/Dropbox/UCR_2020/SP1649-Espacial/0.clases/data/TCE.shp", layer: "TCE"
## with 11 features
## It has 5 fields
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/marce/Dropbox/UCR_2020/SP1649-Espacial/0.clases/data/NY8cities.shp", layer: "NY8cities"
## with 6 features
## It has 1 fields

Datos de casos de leucemia (transformados en escala logarítmica)

\[Z = \frac{1000log(Y_i+1)}{n_i}\]

library(spdep)
nylm <- lm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=NY8)
summary(nylm)
## 
## Call:
## lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7417 -0.3957 -0.0326  0.3353  4.1398 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.51728    0.15856  -3.262  0.00124 ** 
## PEXPOSURE    0.04884    0.03506   1.393  0.16480    
## PCTAGE65P    3.95089    0.60550   6.525 3.22e-10 ***
## PCTOWNHOME  -0.56004    0.17031  -3.288  0.00114 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6571 on 277 degrees of freedom
## Multiple R-squared:  0.1932, Adjusted R-squared:  0.1844 
## F-statistic:  22.1 on 3 and 277 DF,  p-value: 7.306e-13
plot(nylm$fitted.values, nylm$residuals)

¿Puede mejorar, no? Debido a que no podemos asegurar la independencia entre las unidades estadísticas, entonces podemos empezar a construir una estructura espacial que pueda adaptarse a segmentos censales.

¿Cómo?

  1. Vecinos, vecindades y pesos.
  2. Medidas de autocorrelación espacial.
  3. Modelos para medidas dependientes en áreas.
plot(NY8, border="grey60", axes=TRUE)
text(coordinates(cities), labels=as.character(cities$names), font=2, cex=0.9)
text(bbox(NY8)[1,1], bbox(NY8)[2,2], labels="a)", cex=0.8)

plot(NY8, border="grey60", axes=TRUE)
points(TCE, pch=1, cex=0.7)
points(TCE, pch=3, cex=0.7)
text(coordinates(TCE), labels=as.character(TCE$name), cex=0.7,
 font=1, pos=c(4,1,4,1,4,4,4,2,3,4,2), offset=0.3)
text(bbox(NY8)[1,1], bbox(NY8)[2,2], labels="b)", cex=0.8)

plot(NY8, border="grey60", axes=TRUE)

plot(NY8, border="grey60", axes=TRUE)
plot(NY_nb, coordinates(NY8), pch=19, cex=0.6, add=TRUE)

Syracuse <- NY8[NY8$AREANAME == "Syracuse city",]
Sy0_nb <- subset(NY_nb, NY8$AREANAME == "Syracuse city")
summary(Sy0_nb)
## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 346 
## Percentage nonzero weights: 8.717561 
## Average number of links: 5.492063 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 
##  1  1  5  9 14 17  9  6  1 
## 1 least connected region:
## 164 with 1 link
## 1 most connected region:
## 136 with 9 links

Pero en este caso los vecinos ya estaban calculados en un archivo .gal. ¿Cómo podemos crearlos desde el archivo .shp?

coords <- coordinates(Syracuse)
IDs <- row.names(Syracuse)

# Vecinos por movimientos de ajedrez:
Sy1_nb <- poly2nb(Syracuse) # reina default
plot(Sy1_nb, coordinates(Syracuse), pch=19, cex=0.6)

Sy2_nb <- poly2nb(Syracuse, queen=FALSE) # torre
plot(Sy2_nb, coordinates(Syracuse), pch=19, cex=0.6)

## Por distancia:
Sy8_nb <- knn2nb(knearneigh(coords, k=1), row.names=IDs)
plot(Sy8_nb, coordinates(Syracuse), pch=19, cex=0.6)

Sy9_nb <- knn2nb(knearneigh(coords, k=2), row.names=IDs)
plot(Sy9_nb, coordinates(Syracuse), pch=19, cex=0.6)

Sy10_nb <- knn2nb(knearneigh(coords, k=4), row.names=IDs)
plot(Sy10_nb, coordinates(Syracuse), pch=19, cex=0.6)

dsts <- unlist(nbdists(Sy8_nb, coords))
Sy11_nb <- dnearneigh(coords, d1=0, d2=0.75*max(dsts), row.names=IDs)
plot(Sy11_nb, coordinates(Syracuse), pch=19, cex=0.6)

Ya tenemos la definición de nuestros vecinos, según tenga sentido en nuestra pregunta de investigación. Ahora, ¿cómo definir cuál vecino es más o menos importante?

¿Qué es una matriz de pesos W?

https://msarrias.weebly.com/uploads/3/7/7/8/37783629/lecture1.pdf

Estilo W: estandarización por fila. Cada fila debe sumar 1.

Sy0_lw_W <- nb2listw(Sy0_nb)
Sy0_lw_W
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 346 
## Percentage nonzero weights: 8.717561 
## Average number of links: 5.492063 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1      S2
## W 63 3969 63 24.78291 258.564
names(Sy0_lw_W)
## [1] "style"      "neighbours" "weights"
names(attributes(Sy0_lw_W))
## [1] "names"     "class"     "region.id" "call"      "GeoDa"

Revisemos estas sumas:

1/rev(range(card(Sy0_lw_W$neighbours)))
## [1] 0.1111111 1.0000000
summary(unlist(Sy0_lw_W$weights))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1111  0.1429  0.1667  0.1821  0.2000  1.0000
summary(sapply(Sy0_lw_W$weights, sum))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

Estilo B: binario. Unos cuando son vecinos:

Sy0_lw_B <- nb2listw(Sy0_nb, style="B")
summary(unlist(Sy0_lw_B$weights))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1
summary(sapply(Sy0_lw_B$weights, sum))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   4.500   6.000   5.492   6.500   9.000

Otros estilos: C, S y U ver ?nb2listw

“B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme proposed by Tiefelsdorf et al. 1999, p. 167-168 (sums over all links to n).”

La escogencia de criterios de pesos afecta el resultado, por lo que Uds deben tener o: a) una muy buena razón en el contexto de aplicación para usarla o b) evidencia teórica o empírica de que el método funciona mejor para la naturaleza de la aplicación. La última opción (S) tiene evidencia de minimizar ese impacto.

Ejemplo: suponga que creemos que la relación es más fuerte si los vecinos son más cercanos (similar a geoestadística). Entonces, podríamos usar el IDW que ya conocemos. Sin embargo, podríamos considerar distintas relaciones según la dirección (aunque esto complica un supuesto de simetría que veremos más adelante).

dsts <- nbdists(Sy0_nb, coordinates(Syracuse))
idw <- lapply(dsts, function(x) 1/(x/1000))
Sy0_lw_idwB <- nb2listw(Sy0_nb, glist=idw, style="B")
summary(unlist(Sy0_lw_idwB$weights))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3886  0.7374  0.9259  0.9963  1.1910  2.5274
summary(sapply(Sy0_lw_idwB$weights, sum))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.304   3.986   5.869   5.471   6.737   9.435

Ahora que tenemos los vecinos y los pesos, podemos calcular distancias para cada estilo:

library(RColorBrewer)
pal <- brewer.pal(9, "Reds")
oopar <- par(mfrow=c(1,3), mar=c(1,1,3,1)+0.1)
z <- t(listw2mat(Sy0_lw_W))
brks <- c(0,0.1,0.143,0.167,0.2,0.5,1)
nbr3 <- length(brks)-3
image(1:63, 1:63, z[,ncol(z):1], breaks=brks, col=pal[c(1,(9-nbr3):9)],
 main="W style", axes=FALSE)
box()
z <- t(listw2mat(Sy0_lw_B))
image(1:63, 1:63, z[,ncol(z):1], col=pal[c(1,9)], main="B style", axes=FALSE)
box()
z <- t(listw2mat(Sy0_lw_idwB))
brks <- c(0,0.35,0.73,0.93,1.2,2.6)
nbr3 <- length(brks)-3
image(1:63, 1:63, z[,ncol(z):1], breaks=brks, col=pal[c(1,(9-nbr3):9)],
 main="IDW B style", axes=FALSE)
box()

par(oopar)
  1. ¿En cuáles tenemos simetría?
  2. ¿Qué significa el color? Oscuro: cercano a uno, claro cercano a cero.
  3. ¿Qué pasa con los valores de 0? ¿son NAs?
try(Sy0_lw_D1 <- nb2listw(Sy11_nb, style="B"))
## Error in nb2listw(Sy11_nb, style = "B") : Empty neighbour sets found
Sy0_lw_D1 <- nb2listw(Sy11_nb, style="B", zero.policy=TRUE)
print(Sy0_lw_D1, zero.policy=TRUE)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 230 
## Percentage nonzero weights: 5.794911 
## Average number of links: 3.650794 
## 2 regions with no links:
## 154 168
## 
## Weights style: B 
## Weights constants summary:
##    n   nn  S0  S1   S2
## B 61 3721 230 460 4496

El problema de datos sin vecinos se puede resolver con complete.cases, o subset.listw.

Una vez creados los vecindarios, pesos y distancias, ¿cómo las guardamos?

Sy14_nb <- read.gal("data/Sy_GeoDa1.GAL")
isTRUE(all.equal(Sy0_nb, Sy14_nb, check.attributes=FALSE))
## [1] TRUE
Sy16_nb <- read.gwt2nb("data/Sy_GeoDa4.GWT")
## Warning in read.gwt2nb("data/Sy_GeoDa4.GWT"): Old-style GWT file
isTRUE(all.equal(Sy10_nb, Sy16_nb, check.attributes=FALSE))
## [1] TRUE
## esto es muy útil para Winbugs, Stan, INLA:
library(foreign)
df <- as.data.frame(listw2mat(Sy0_lw_B))
write.dta(df, file="Sy0_lw_B.dta", version=7)

Usando los pesos para simular estructuras de autocorrelación espacial:

\((I - \rho W)\) usando la función invIrW:

set.seed(987654)
n <- length(Sy0_nb)
uncorr_x <- rnorm(n)
rho <- 0.5
autocorr_x <- invIrW(Sy0_lw_W, rho) %*% uncorr_x
## Registered S3 methods overwritten by 'spatialreg':
##   method                   from 
##   residuals.stsls          spdep
##   deviance.stsls           spdep
##   coef.stsls               spdep
##   print.stsls              spdep
##   summary.stsls            spdep
##   print.summary.stsls      spdep
##   residuals.gmsar          spdep
##   deviance.gmsar           spdep
##   coef.gmsar               spdep
##   fitted.gmsar             spdep
##   print.gmsar              spdep
##   summary.gmsar            spdep
##   print.summary.gmsar      spdep
##   print.lagmess            spdep
##   summary.lagmess          spdep
##   print.summary.lagmess    spdep
##   residuals.lagmess        spdep
##   deviance.lagmess         spdep
##   coef.lagmess             spdep
##   fitted.lagmess           spdep
##   logLik.lagmess           spdep
##   fitted.SFResult          spdep
##   print.SFResult           spdep
##   fitted.ME_res            spdep
##   print.ME_res             spdep
##   print.lagImpact          spdep
##   plot.lagImpact           spdep
##   summary.lagImpact        spdep
##   HPDinterval.lagImpact    spdep
##   print.summary.lagImpact  spdep
##   print.sarlm              spdep
##   summary.sarlm            spdep
##   residuals.sarlm          spdep
##   deviance.sarlm           spdep
##   coef.sarlm               spdep
##   vcov.sarlm               spdep
##   fitted.sarlm             spdep
##   logLik.sarlm             spdep
##   anova.sarlm              spdep
##   predict.sarlm            spdep
##   print.summary.sarlm      spdep
##   print.sarlm.pred         spdep
##   as.data.frame.sarlm.pred spdep
##   residuals.spautolm       spdep
##   deviance.spautolm        spdep
##   coef.spautolm            spdep
##   fitted.spautolm          spdep
##   print.spautolm           spdep
##   summary.spautolm         spdep
##   logLik.spautolm          spdep
##   print.summary.spautolm   spdep
##   print.WXImpact           spdep
##   summary.WXImpact         spdep
##   print.summary.WXImpact   spdep
##   predict.SLX              spdep

Los datos con y sin correlación y con autocorrelación se ven así:

oopar <- par(mfrow=c(1,2), mar=c(4,4,3,2)+0.1)
plot(uncorr_x, lag(Sy0_lw_W, uncorr_x), xlab="", cex.lab=0.8,
 ylab="spatial lag", main="Uncorrelated random variable", cex.main=0.8)
lines(lowess(uncorr_x, lag(Sy0_lw_W, uncorr_x)), lty=2, lwd=2)
plot(autocorr_x, lag(Sy0_lw_W, autocorr_x),
 xlab="", ylab="",
 main="Autocorrelated random variable", cex.main=0.8, cex.lab=0.8)
lines(lowess(autocorr_x, lag(Sy0_lw_W, autocorr_x)), lty=2, lwd=2)

par(oopar)

Pruebas para autocorrelación espacial:

Moran I test: Antes de hacer cualquier tipo de test de autocorrelación espacial, debemos asegurarnos (de la misma manera que con series de tiempo) que no existe ningún patrón espacial en la media. Para ello, puede modelarse la media utilizando otras variables, y luego hacer el test con los residuales. También, supuestos acerca de la distribución de los datos deben estar presentes en este caso (hemos visto como asumimos casi siempre normalidad, por ejemplo).

Fórmula:

\(I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}}\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(y_i-\bar{y})(y_j-\bar{y})}{\sum_{i=1}^{n}\sum_{j=1}^{n}(y_j-\bar{y})^2}\)

donde \(y_i\) es la observación i, \(\bar{y}\) es la media de la variable de interés y \(w_{ij}\) es el peso definido entre la observación i y la j. Utilizar la media en este caso, implica que la media es la misma para todas las observaciones en el espacio, es decir, que no existen patrones.

Ahora vamos a calcular el test para varios casos:

moran.test(uncorr_x, listw=Sy0_lw_W)
## 
##  Moran I test under randomisation
## 
## data:  uncorr_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = -0.22698, p-value = 0.5898
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.033287124      -0.016129032       0.005714128
moran.test(autocorr_x, listw=Sy0_lw_W)
## 
##  Moran I test under randomisation
## 
## data:  autocorr_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = 3.1027, p-value = 0.0009588
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.218178182      -0.016129032       0.005702793
moran.test(autocorr_x, listw=nb2listw(Sy9_nb, style="W"))
## 
##  Moran I test under randomisation
## 
## data:  autocorr_x  
## weights: nb2listw(Sy9_nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.8619, p-value = 0.03131
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.19207401       -0.01612903        0.01250486
et <- coords[,1] - min(coords[,1])
trend_x <- uncorr_x + 0.00025 * et
moran.test(trend_x, listw=Sy0_lw_W)
## 
##  Moran I test under randomisation
## 
## data:  trend_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = 3.3438, p-value = 0.0004131
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.237474257      -0.016129032       0.005751989
lm.morantest(lm(trend_x ~ et), listw=Sy0_lw_W)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = trend_x ~ et)
## weights: Sy0_lw_W
## 
## Moran I statistic standard deviate = -0.31166, p-value = 0.6223
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     -0.053830651     -0.030937453      0.005395814

Resultados:

  1. Moran’s I: Valor entre -1 y 1.
  2. \(E(I) = \frac{-1}{n-1}\) para los casos centrados en la media.
  3. La variancia se aproxima por medio de remuestreo.
  4. \(St. deviate = \frac{I-E(I)}{\sqrt{var(I)}}\)

Los valores significativamente por debajo de -1 / (N-1) indican una autocorrelación espacial negativa y los valores por encima de -1 / (N-1) indican una autocorrelación espacial positiva. Para las pruebas de hipótesis estadísticas, los valores I de Moran pueden transformarse en puntajes z.

Ojo con las conclusiones, ¿cómo cambian con los mismos datos pero definiendo vecinos de otra manera o modelando la media distinto?

Test globales:

Existen otros tipos de tests como Geary’s C, Getis-Ord G, Mantest test, etc. Sin embargo, el test I de Moran, está aceptado como referencia en el caso de estadística de áreas. Todos son para variables continuas, pero existen correcciones para variables categóricas (para el test de Moran).

Veamos otros ejemplos usando Moran. Mismos datos, pero usando una estructura de vecinos distinta:

K <- moran(NY8$Cases, listw=nb2listw(NY_nb, style="B"), n=length(NY8$Cases), S0=Szero(nb2listw(NY_nb, style="B")))$K
moran.test(NY8$Cases, listw=nb2listw(NY_nb)) # default es W
## 
##  Moran I test under randomisation
## 
## data:  NY8$Cases  
## weights: nb2listw(NY_nb)    
## 
## Moran I statistic standard deviate = 3.9778, p-value = 3.477e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.146882990      -0.003571429       0.001430595
lw_B <- nb2listw(NY_nb, style="B")
moran.test(NY8$Cases, listw=lw_B)
## 
##  Moran I test under randomisation
## 
## data:  NY8$Cases  
## weights: lw_B    
## 
## Moran I statistic standard deviate = 3.1862, p-value = 0.0007207
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.110387402      -0.003571429       0.001279217

Ojo al cambio en el valor p. 

También, podemos asumir normalidad para calcular la variancia:

moran.test(NY8$Cases, listw=lw_B, randomisation=FALSE)
## 
##  Moran I test under normality
## 
## data:  NY8$Cases  
## weights: lw_B    
## 
## Moran I statistic standard deviate = 3.1825, p-value = 0.0007301
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.110387402      -0.003571429       0.001282228
lm.morantest(lm(Cases ~ 1, NY8), listw=lw_B)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Moran I statistic standard deviate = 3.1825, p-value = 0.0007301
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.110387402     -0.003571429      0.001282228

Detalles: el test de Moran para bajo el supuesto de normalidad es el mismo que el test de residuos (lm(formula = Cases ~ 1, data = NY8)).

También, podemos usar la versión exacta o la aproximación de Saddlepoint (punto de silla). El primero asume normalidad y el segundo no. En este caso no hay mucha diferencia, pero si la puede haber cuando no se cumple la normalidad.

lm.morantest.sad(lm(Cases ~ 1, NY8), listw=lw_B)
## 
##  Saddlepoint approximation for global Moran's I (Barndorff-Nielsen
##  formula)
## 
## data:  
## model:lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Saddlepoint approximation = 2.9929, p-value = 0.001382
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I 
##        0.1103874
lm.morantest.exact(lm(Cases ~ 1, NY8), listw=lw_B)
## 
##  Global Moran I statistic with exact p-value
## 
## data:  
## model:lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Exact standard deviate = 2.9923, p-value = 0.001384
## alternative hypothesis: greater
## sample estimates:
## [1] 0.1103874

Otra manera de calcular la I de Moran es utilizando un test de permutación: ¿qué pasa si asigno locaciones aleatorias a valores? ¿obtengo el mismo estadístico o este cambia? La hipótesis nula es que tienen la misma distribución. Si tengo evidencia para rechazar la hipótesis nula, entonces significa que la locación original importa y que por lo tanto si existe autocorrelación espacial.

set.seed(1234)
bperm <- moran.mc(NY8$Cases, listw=lw_B, nsim=999)
bperm
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  NY8$Cases 
## weights: lw_B  
## number of simulations + 1: 1000 
## 
## statistic = 0.11039, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Otra opción es utilizar un boostrap paramétrico en el caso de tener cuentas en lugar de variables aleatorias contínuas:

r <- sum(NY8$Cases)/sum(NY8$POP8) # proporción de casos global
rni <- r*NY8$POP8 # crea un valor esperado para cada caso.
CR <- function(var, mle) rpois(length(var), lambda=mle) # genera poissons
MoranI.pboot <- function(var, i, listw, n, S0, ...) {
  return(moran(x=var, listw=listw, n=n, S0=S0)$I)
}
set.seed(1234)
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
## 
##     melanoma
boot2 <- boot(NY8$Cases, statistic=MoranI.pboot, R=999, sim="parametric",
  ran.gen=CR, listw=lw_B, n=length(NY8$Cases), S0=Szero(lw_B), mle=rni)
pnorm((boot2$t0 - mean(boot2$t))/sd(boot2$t[,1]), lower.tail=FALSE)
## [1] 0.1472112

¿Cómo se ve?

oopar <- par(mfrow=c(1,2))
xlim <- range(c(bperm$res, boot2$t[,1]))
hist(bperm$res[-length(bperm$res)], main="Permutation bootstrap", xlab=expression(I[std]), xlim=xlim, density=15, angle=45, ylim=c(0,260))
abline(v=bperm$statistic, lty=2)
hist(boot2$t, col=rgb(0.4,0.4,0.4), main="Parametric bootstrap", xlab=expression(I[CR]), xlim=xlim, ylim=c(0,260))
hist(bperm$res[-length(bperm$res)], density=15, angle=45, add=TRUE)
abline(v=boot2$t0, lty=2)

par(oopar)
## rni <- fitted(glm(Cases ~ 1 + offset(log(POP8)), data=NY8, family="poisson"))

¿A cuál le hacemos caso?

Cuando tenemos cuentas pero con poblaciones muy desiguales entre un área y otra, es recomendable utilizar el EBI Moran:

set.seed(1234)
EBImoran.mc(n=NY8$Cases, x=NY8$POP8, listw=nb2listw(NY_nb, style="B"), nsim=999)
## The default for subtract_mean_in_numerator set TRUE from February 2016
## 
##  Monte-Carlo simulation of Empirical Bayes Index (mean subtracted)
## 
## data:  cases: NY8$Cases, risk population: NY8$POP8
## weights: nb2listw(NY_nb, style = "B")
## number of simulations + 1: 1000
## 
## statistic = 0.072497, observed rank = 983, p-value = 0.017
## alternative hypothesis: greater

Es decir, una vez que tomamos en cuenta las distintas poblaciones, si podemos detectar autocorrelación.

Esta función puede ayudarnos a detectar donde está la autocorrelación:

cor8 <- sp.correlogram(neighbours=NY_nb, var=NY8$Cases, order=8, method="I", style="C")
library(pgirmess)
corD <- correlog(coordinates(NY8), NY8$Cases, method="Moran")
oopar <- par(mfrow=c(1,2))
plot(cor8, main="Contiguity lag orders")
plot(corD, main="Distance bands")

par(oopar)

Test Locales: Podemos usar las mismas herramientas para detectar clusters o agrupaciones.

\(I_i = \frac{(y_i-\bar{y})\sum_{j=1}^{n}w_{ij}(y_j-\bar{y})}{\frac{\sum_{i=1}^n (y_i-\bar{y})^2}{n}}\)

oopar <- par(mfrow=c(1,2))
msp <- moran.plot(NY8$Cases, listw=nb2listw(NY_nb, style="C"), quiet=TRUE)
title("Moran scatterplot")
infl <- apply(msp["is_inf"], 1, any)
x <- NY8$Cases
lhx <- cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=TRUE)
wx <- lag(nb2listw(NY_nb, style="C"), NY8$Cases)
lhwx <- cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=TRUE)
lhlh <- interaction(lhx, lhwx, infl, drop=TRUE)
cols <- rep(1, length(lhlh))
cols[lhlh == "H.L.TRUE"] <- 2
cols[lhlh == "L.H.TRUE"] <- 3
cols[lhlh == "H.H.TRUE"] <- 4
plot(NY8, col=brewer.pal(4, "Accent")[cols])
legend("topright", legend=c("None", "HL", "LH", "HH"), fill=brewer.pal(4, "Accent"), bty="n", cex=0.8, y.intersp=0.8)
title("Tracts with influence")

par(oopar)

Esto no es una tarea fácil, por lo que más pruebas están a nuestra disposición:

## moran.plot(NY8$Cases, listw=nb2listw(NY_nb, style="C"))
lm1 <- localmoran(NY8$Cases, listw=nb2listw(NY_nb, style="C"))
lm2 <- as.data.frame(localmoran.sad(lm(Cases ~ 1, NY8), nb=NY_nb, style="C"))
lm3 <- as.data.frame(localmoran.exact(lm(Cases ~ 1, NY8), nb=NY_nb, style="C"))
r <- sum(NY8$Cases)/sum(NY8$POP8)
rni <- r*NY8$POP8
lw <- nb2listw(NY_nb, style="C")
sdCR <- (NY8$Cases - rni)/sqrt(rni)
wsdCR <- lag(lw, sdCR)
I_CR <- sdCR * wsdCR

Ejemplo: cálculo del Test de Moran local para cada área, bajo el supuesto standard, y bajo el supuesto de riesgo constante (corrección por distinta población en cada área):

library(RColorBrewer)
gry <- c(rev(brewer.pal(8, "Reds")[1:6]), brewer.pal(6, "Blues"))
NY8$Standard <- lm1[,1]
NY8$"Constant_risk" <- I_CR
#nms <- match(c("Standard", "Constant_risk"), names(NY8))
spplot(NY8, c("Standard", "Constant_risk"), at=c(-2.5,-1.4,-0.6,-0.2,0,0.2,0.6,4,7), col.regions=colorRampPalette(gry)(8))

Podríamos también encontrar la probabilidad (por medio de un test de Monte Carlo) para cada supuesto. Valores cercanos a cero indican agrupaciones en los datos (clusters: valores de áreas cercanas se parecen) y valores cercanos a uno indican hotspots: valores en un área difieren mucho comparadas con sus vecinos.s

set.seed(1234)
nsim <- 999
N <- length(rni)
sims <- matrix(0, ncol=nsim, nrow=N)
for (i in 1:nsim) {
  y <- rpois(N, lambda=rni)
  sdCRi <- (y - rni)/sqrt(rni)
  wsdCRi <- lag(lw, sdCRi)
  sims[,i] <- sdCRi * wsdCRi 
}
xrank <- apply(cbind(I_CR, sims), 1, function(x) rank(x)[1])
diff <- nsim - xrank
diff <- ifelse(diff > 0, diff, 0)
pval <- punif((diff + 1)/(nsim + 1))

NY8$Normal <- lm2[,3]
NY8$Randomisation <- lm1[,5]
NY8$Saddlepoint <- lm2[,5]
NY8$Exact <- lm3[,5]
NY8$Constant_risk <- pval
gry <- c(rev(brewer.pal(6, "Reds")), brewer.pal(6, "Blues"))
spplot(NY8, c("Normal", "Randomisation", "Saddlepoint", "Exact", "Constant_risk"), at=c(0,0.01,0.05,0.1,0.9,0.95,0.99,1), col.regions=colorRampPalette(gry)(7))

spplot(NY8, c("Normal", "Exact", "Constant_risk"), xlim=c(405200, 432200), ylim=c(4652700, 4672000), at=c(0,0.01,0.05,0.1,0.9,0.95,0.99,1), col.regions=colorRampPalette(gry)(7))

oopar <- par(mfrow=c(1,2))
msp <- moran.plot(NY8$Cases, listw=nb2listw(NY_nb, style="W"), quiet=TRUE)
title("Moran scatterplot")
infl <- apply(msp["is_inf"], 1, any)
x <- NY8$Cases
lhx <- cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=TRUE)
wx <- lag(nb2listw(NY_nb, style="C"), NY8$Cases)
lhwx <- cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=TRUE)
lhlh <- interaction(lhx, lhwx, infl, drop=TRUE)
cols <- rep(1, length(lhlh))
cols[lhlh == "H.L.TRUE"] <- 2
cols[lhlh == "L.H.TRUE"] <- 3
cols[lhlh == "H.H.TRUE"] <- 4
plot(NY8, col=brewer.pal(4, "Accent")[cols])
legend("topright", legend=c("None", "HL", "LH", "HH"), fill=brewer.pal(4, "Accent"), bty="n", cex=0.8, y.intersp=0.8)
title("Tracts with influence")

Laboratorio 5:

Preguntas 1 al 5 del cuaderno de trabajo, con una breve introducción acerca del origen de los datos, y el objetivo del análisis: https://rspatial.org/raster/analysis/3-spauto.html